R Notebook

---

Začneme tým, že si do RStudia importujeme dáta. Ide o surové textové dáta, takže to nie je žiadny problém. Zvolíme import dát, zadáme cestu a prenesieme do nášho kódu:

hygdata_v3 <- read.csv("~/GitHub/HYG-Database/hygdata_v3.csv")
View(hygdata_v3)

Teraz už s nimi môžeme pracovať, najprv sa však pozrieme, ako sú usporiadané.

head(hygdata_v3)

Môžeme si všimnúť, že ide o pomerne komplexné dáta, ktoré sú bohaté okrem zaradenia do rôznych astronomických katalógov aj na údaje o relatívnej polohe a rýchlosti týchto vesmírnych objektov, ako aj o ich svietivosti a podobne. Máme teda na výber množstvo súvislostí, ktoré môžeme odhaľovať medzi jednotlivými parametrami.

Našou úplne prvou úlohou bude pozrieť sa, ako sú naše dáta v súbore rozdelené. Nebudeme sa trápiť s vykresľovaním jednotlivých charakteristík, radšej použijeme nástroje, ktoré sú v R už vbudované (summary). Zároveň by bolo ale fajn otestovať nedokonalosti v dátach, a ďalšiu informáciu nám môže poskytnúť aj modus, ktorý v bežnom summary nie je prítomný, a tak si ho pripravíme zvlášť a aplikujeme spolu so summary na naše dáta.

print(summary(hygdata_v3))

       id              hip               hd               hr                gl                  bf        
 Min.   :     0   Min.   :     1   Min.   :     1   Min.   :   2             :115813             :116516  
 1st Qu.: 29903   1st Qu.: 29565   1st Qu.: 46734   1st Qu.:2281     GJ 1001 :     1   41The1Ori :     3  
 Median : 59807   Median : 59172   Median :110382   Median :4567     GJ 1002 :     1   100    Her:     2  
 Mean   : 59807   Mean   : 59170   Mean   :114377   Mean   :4564     GJ 1003 :     1   12Omi Cap :     2  
 3rd Qu.: 89710   3rd Qu.: 88763   3rd Qu.:175849   3rd Qu.:6850     GJ 1004 :     1   17Sig CrB :     2  
 Max.   :119615   Max.   :120404   Max.   :358431   Max.   :9110     GJ 1005A:     1   19    Lyn :     2  
                  NA's   :1658     NA's   :20732    NA's   :110589   (Other) :  3796   (Other)   :  3087  
        proper             ra              dec               dist               pmra          
           :119468   Min.   : 0.000   Min.   :-89.782   Min.   :     0.0   Min.   :-4432.650  
 268 G. Cet:     1   1st Qu.: 6.217   1st Qu.:-36.421   1st Qu.:   115.1   1st Qu.:  -15.460  
 3C 273    :     1   Median :12.128   Median : -1.640   Median :   213.7   Median :   -1.680  
 82 G. Eri :     1   Mean   :12.095   Mean   : -1.985   Mean   :  8764.8   Mean   :   -1.312  
 96 G. Psc :     1   3rd Qu.:18.117   3rd Qu.: 31.519   3rd Qu.:   390.6   3rd Qu.:   12.180  
 Acamar    :     1   Max.   :23.999   Max.   : 89.569   Max.   :100000.0   Max.   : 6767.260  
 (Other)   :   141                                                                            
     pmdec                rv                 mag             absmag             spect      
 Min.   :-5813.00   Min.   :-386.9000   Min.   :-26.70   Min.   :-16.6800   K0     : 8569  
 1st Qu.:  -22.39   1st Qu.:   0.0000   1st Qu.:  7.65   1st Qu.:  0.1390   G5     : 6021  
 Median :   -5.76   Median :   0.0000   Median :  8.46   Median :  1.4950   A0     : 4943  
 Mean   :  -19.32   Mean   :  -0.2773   Mean   :  8.43   Mean   :  0.9921   F8     : 4366  
 3rd Qu.:    3.77   3rd Qu.:   0.0000   3rd Qu.:  9.17   3rd Qu.:  3.1590   G0     : 4237  
 Max.   : 9999.99   Max.   : 471.0000   Max.   : 21.00   Max.   : 19.6290   F5     : 3860  
                                                                            (Other):87618  
       ci                x                   y                   z                   vx            
 Min.   :-0.4000   Min.   :-99950.39   Min.   :-99979.25   Min.   :-99964.98   Min.   :-1.023e-01  
 1st Qu.: 0.3480   1st Qu.:   -89.02   1st Qu.:   -91.17   1st Qu.:  -107.54   1st Qu.:-1.034e-05  
 Median : 0.6160   Median :    -1.04   Median :    -1.24   Median :    -3.41   Median : 1.300e-07  
 Mean   : 0.7115   Mean   :  -234.85   Mean   :   -39.57   Mean   :  -231.90   Mean   :-2.891e-05  
 3rd Qu.: 1.0830   3rd Qu.:    86.27   3rd Qu.:    91.85   3rd Qu.:    94.97   3rd Qu.: 9.100e-06  
 Max.   : 5.4600   Max.   : 99982.37   Max.   : 99996.07   Max.   : 99862.51   Max.   : 1.023e-01  
 NA's   :1882                                                                                      
       vy                   vz                 rarad           decrad            pmrarad          
 Min.   :-1.023e-01   Min.   :-1.023e-01   Min.   :0.000   Min.   :-1.56700   Min.   :-2.149e-05  
 1st Qu.:-1.860e-06   1st Qu.:-1.998e-05   1st Qu.:1.628   1st Qu.:-0.63566   1st Qu.:-7.495e-08  
 Median : 1.182e-05   Median :-6.230e-06   Median :3.175   Median :-0.02862   Median :-8.140e-09  
 Mean   : 2.164e-04   Mean   :-1.566e-04   Mean   :3.167   Mean   :-0.03465   Mean   :-6.360e-09  
 3rd Qu.: 3.209e-05   3rd Qu.: 3.150e-06   3rd Qu.:4.743   3rd Qu.: 0.55011   3rd Qu.: 5.905e-08  
 Max.   : 1.023e-01   Max.   : 1.023e-01   Max.   :6.283   Max.   : 1.56328   Max.   : 3.281e-05  
                                                                                                  
    pmdecrad              bayer             flam             con             comp        comp_primary   
 Min.   :-2.818e-05          :118078   Min.   :  1.0    Cen    : 4222   Min.   :1.000   Min.   :     0  
 1st Qu.:-1.086e-07   Alp    :    80   1st Qu.: 14.0    UMa    : 3549   1st Qu.:1.000   1st Qu.: 29812  
 Median :-2.793e-08   Bet    :    77   Median : 31.5    Her    : 3383   Median :1.000   Median : 59628  
 Mean   :-9.366e-08   Eps    :    74   Mean   : 37.3    Cyg    : 3081   Mean   :1.005   Mean   : 59635  
 3rd Qu.: 1.828e-08   Del    :    71   3rd Qu.: 55.0    Hya    : 3027   3rd Qu.:1.000   3rd Qu.: 89453  
 Max.   : 5.007e-05   Eta    :    68   Max.   :140.0    Cet    : 2958   Max.   :3.000   Max.   :119615  
                      (Other):  1166   NA's   :116874   (Other):99394                                   
       base             lum                 var            var_min          var_max      
         :118527   Min.   :        0          :113624   Min.   :-1.33    Min.   :-1.52   
 GJ 1047 :     3   1st Qu.:        5   R      :    60   1st Qu.: 8.53    1st Qu.: 8.24   
 Gl  57.1:     3   Median :       22   S      :    52   Median : 9.85    Median : 9.65   
 Gl  60  :     3   Mean   :   354469   T      :    50   Mean   : 9.50    Mean   : 9.26   
 Gl 100  :     3   3rd Qu.:       77   W      :    39   3rd Qu.:10.71    3rd Qu.:10.49   
 Gl 106.1:     3   Max.   :409260660   U      :    37   Max.   :14.90    Max.   :13.70   
 (Other) :  1072                       (Other):  5752   NA's   :102623   NA's   :102623  

Ak si pozorne prezrieme zobrazené charakteristiky, môžeme si všimnúť, že v niektorých položkách sú prítomné outlieri. Keďže však ide o reálne objekty s objektívnymi vlastnosťami, nebudeme ich zatiaľ vyradzovať (napríklad lum - svietivosť ako násobok svietivosti Slnka). Určite však budeme musieť upraviť veľmi dôležitú položku dist - teda vzdialenosť hviezdy získanú na základe jej paralaxy (relatívnej zmeny polohy kvôli pozorovaniu z dvoch protiľahlých miest na obežnej dráhe Slnka). Ako si môžeme všimnúť, tretí kvartil zasahuje hlboko pod maximálnu hodnotu súboru. Keď sa pozrieme na modus tohto spojitého súboru, zistíme, že je to práve táto maximálna hodnota. Táto maximálna presná hodnota v spojitom spektre je dosť podozrivá:

modus <- function(v) 
{
   uniqv <- unique(v)
   uniqv[which.max(tabulate(match(v, uniqv)))]
}
modus(hygdata_v3[,10])

[1] 1e+05

hist(hygdata_v3[,10],xlab = "vzdialenosť(získaná pomocou paralaxy)",ylab = "počet výskytov",col = "darkblue",border = "blue")

Aha naozaj. Nájsť modus v spojitom súbore dát väčšinou znamená nejakú špecifickú vlastnosť značenia a histogram to potvrdzuje. V tomto bode si teda ešte raz prezrieme readme ku našej databáze. Naozaj. hodnota 100 000 znamená, že dáta chýbajú, alebo sú pochybné. Očistime naše dáta od tohoto znečistenia a pozrime sa na početnosti.

library(tidyverse)

-- Attaching packages --------------------------------------- tidyverse 1.2.1 --
<U+221A> ggplot2 3.2.1     <U+221A> purrr   0.3.3
<U+221A> tibble  2.1.3     <U+221A> dplyr   0.8.3
<U+221A> tidyr   1.0.0     <U+221A> stringr 1.4.0
<U+221A> readr   1.3.1     <U+221A> forcats 0.4.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

purified_1 <- filter(hygdata_v3,dist<=50000)
rest_1 <- filter(hygdata_v3, dist > 50000)
hist(purified_1[,10],xlab = "vzdialenosť(získaná pomocou paralaxy)",ylab = "počet výskytov",col = "darkblue",border = "blue")

hist(rest_1[,14],xlab = "magnitúda",ylab = "počet výskytov",col = "darkred",border = "gray")

Začnime s niečím jednoduchým. Radi by sme dokázali, že magnitúda (zdanlivá jasnosť hviezd) bude rásť (jej reprezentujúca číslica) so vzdialenosťou (prirodzený predpoklad v prípade homogénneho prostredia). Vynesme si teda tieto dve veličiny do grafu.

library(ggplot2)
ggplot(data = purified_1) +
geom_point(mapping = aes(x = dist, y = mag,size = 1/10^(mag),alpha = 1/100))

Z takéhoto grafu až tak veľa informácií nezískame. Dáta by bolo fajn si utriediť do nejakých skupín, a až takéto skupiny zobrazovať (a prípadne nastaviť vizuál podľa ich veľkostí). Poďme sa teda pustiť do našej prvej podúlohy - vyberme si z dát len zaujímavé položky a zoraďme ich podľa hodnoty.

hygdata_1 <- select(purified_1, id, proper, dist, rv, mag, absmag, con, lum)
hygdata_1

Teraz sa pokúsime nejako regulovať počet údajov. Jednotlivé objekty podľa vzdialenosti a magnitúdy zatriedime do vhodného počtu skupín. Využijeme zoradenie objektov vzostupne podľa vzdialenosťi, vykreslenie príslušnosti do skupín a funkciu cut.

dist

kat <- cut(hygdata_1$dist, 40, include.lowest=TRUE)
mutate(hygdata_1, kat)

hygdata_1

ggplot(data = hygdata_1) +
geom_point(mapping = aes(x = kat, y = mag, size = 1/10^(mag), alpha = 1/100), )

Ako si môžeme všimnúť, stále informačná hodnota takejto vizualizácie, nespĺňa naše požiadavky. Skúsime to inak: vizualizujme priamo závislosť vzdialenosti od magnitúry, ale príslušnosť ku jednotlivým vzialenostiam pre lepšiu prehľadnosť ofarbíme podľa rozdelenia spektra vzdialenosté do nastaviteľného počtu intervalov.

hygdata_1

kat <- cut(hygdata_1$dist, 25, include.lowest=TRUE, labels = FALSE)
mutate(hygdata_1, kat)

ggplot(data = hygdata_1) +
geom_point(mapping = aes(x = dist, y = mag, size = 1/10^(mag), alpha = 1/100), colour = kat)

Ako si môžeme všimnúť, vo funkcií cut sme použili labels = FALSE, pretože nechceme ručne zadávať potenciálne veľké množstvo názvov, ktoré nemajú informačnú hodnotu, len sú feed-om pre parameter colour pri kreslení grafu. To, že sme túto možnosť dali do kódu nám ale umožňuje voľne a nezáväzne adjustovať počet farebných hladín a počet “chlievikov”, ktoré zobrazujeme - vlastne to používame ako taký “kvázihistogram” - vidíme takto na našom grafe početnosti. Táto funkcia to delí na rovnaké podintervaly. Nám by pre vizualizáciu možno pomohlo, keby to delíme na kúsky s rovnakým počtom vzoriek - tak by sme získali informáciu o rozdelení našich dát a bolo by toto rozdelenie krásne vidno. Teraz sa to pokúsime spraviť. Hľadali sme možnosti až sme došli k jednoduchej variante príkazu cut, a to cut_number:

hygdata_1

kat <- cut_number(hygdata_1$dist, 20 , labels = FALSE)
mutate(hygdata_1, kat)

ggplot(data = hygdata_1) +
geom_point(mapping = aes(x = dist, y = mag, size = 1/1.1^(mag), alpha = 1/1000), colour = kat)

Teraz už sú dáta zobrazené pomerne prehľadne. Môžeme teda usúdiť, že náš predpoklad sa nepotvrdil. Vecou, ktorá je vďaka tejto vizualizácií najzrejmejšia, je rozdiel počtov hviezd v databáze pre jednotlivé vzdialenosti. Šírka stĺpcu (chlievika), ktorý obsahuje fixný počet hviezd sa so stúpajúcou vzdialenosťou zväčšuje - nie sme schopní detekovať všetky hviezdy (za predpokladu relatívnej homogenity rozmiestenenia hviezd vo vesmíre). Druhým viditeľným poznatkom je, že do vzdialenosti asi 50 parsekov sme schopní zachytiť aj objekty s veľmi malou zdanlivou magnitúdou. Potom ale počet takýchto objektov v databáze exponenciálne klesá. Tento pokles nie je spôsobený tým, že by vo vzdialenejších častiach menej jasné objekty neboli, ale tým, že ich nie sme schopní detekovať. V našej databáze teda tieto objekty nie sú. Teda ako to už býva, náš jednoduchý predpoklad “so stúpajúcou vzdialenosťou stúpa číslo magnitúdy” musíme korigovať ako “so stúpajúcou vzdialenosťou sme schopní zachytiť čoraz menej hviezd a žiadne s príliš veľkou magnitúdou (nízkou zdanlivou jasnosťou)”.

Druhým cieľom, ktorý by sme chceli splniť, je vykreslenie tvaru súhvezdí na oblohe podľa ich súradníc. Na takýto komplexnejšií problém využijeme fazety a identifikáciu podľa stĺpca con (skratka pre constellation). V prvom rade si z nášho komplexného dátového súboru vyberieme len položky, ktoré sú pre nás zaujímavé. Budeme postupovať analogicky ako v prvej časti.

hygdata_2 <- select(hygdata_v3, id, proper, mag, con, x, y, z, vx, vy, vz, rarad, decrad, pmrarad, pmdecrad)
hygdata_2 <- filter(hygdata_2,con!="")

Poďme sa pozrieť na rozdelenia jednotlivých súradníc - na základe tohoto sa rozhodneme, ktorú metódu vizualizácie zvolíme. Zároveň sa poďme pozrieť, koľko hodnôt môže nadobúdať kategorická premenná con - koľko je súhvezdí?

library(ggplot2)
ggplot(data.frame(hygdata_2$con), aes(x=hygdata_2$con), colour = hygdata_2$con ) +
  geom_bar()

Ako môžeme vidieť, Súhvezdí je trochu moc, aby sme ich mohli všetky vykresliť. Vyberieme si teda len niektoré z nich a tejto vzorke sa budeme venovať. Samozrejme, v prípade záujmu by sme nasledujúce postupy mohli použiť aj na iných ľubovoľných súhvezdiach. Poďme sa pozrieť na počty hviezd v niektorých súhvezdiach / zvoľme napríklad Pisces - Ryby a Leo - lev.

hygdata_2Psc<-filter(hygdata_2, con =="Psc") 
hist(hygdata_2Psc[,3],xlab = "magnitúda",ylab = "počet výskytov",col = "darkred",border = "gray")

hygdata_2Leo<-filter(hygdata_2, con=="Leo") 
hist(hygdata_2Leo[,3],xlab = "magnitúda",ylab = "počet výskytov",col = "darkblue",border = "gray")

Ako vidíme, počty hviezd v jednotlivých súhvezdiach sú vyššie, ako by bolo možné prehľadne zobraziť. Poďme teda odfiltrovať hviezdy, ktoré nie sú viditeľné voľným okom. Uvádza sa, že voľným okom je možné vidieť objekty s veľkosťou okolo 6 magnitúdy (http://dust.ess.uci.edu/ppr/ppr_Wea47.pdf). Použime tento predpoklad a odfiltrujme hviezdy, ktoré nie sú viditeľné voľným okom (po optimalizácií sme zvolili ešte nižšiu hodnotu magnitúdy cutoff):

hygdata_2

hygdata_2_vi <- filter(hygdata_2, mag<4.5)
hygdata_2_vi

A môžeme začať vykresľovať! Samozrejme, ešte zvoľme, ktoré súhvezdia vykreslíme. Tak spravme túto úlohu pre 12 znamení zverokruhu:

hygdata_2_vi

hygdata_2_vi_z <- filter(hygdata_2_vi, con =="Ari"|con =="Tau"|con =="Gem"|con =="Cnc"|con =="Leo"|con =="Vir"|con =="Lib"|con =="Sco"|con =="Sgr"|con =="Aqr"|con =="Psc"|con =="Cap")
hygdata_2_vi_z

A teraz už naozaj môžeme začať vykresľovať. Skúsme, ako to vôbec bude vyzerať pre jedno súhvezdie (dobre identifikovateľný Leo). Vykresľujeme jednotlivé súhvezdia podľa ich súradníc na nebeskej klenbe (najprv sme experimentovali s x, y a z súradnicami až pokým sme si neuvedomili ich význam v tejto databáze). Takto sú charakteristické tvary viditeľné:

hygdata_2_vi_leo <- filter(hygdata_2_vi_z, con =="Leo")
ggplot(data = hygdata_2_vi_leo) +
  geom_point(mapping = aes(x = hygdata_2_vi_leo$rarad, y = hygdata_2_vi_leo$decrad, size = 1/1.01^(hygdata_2_vi_leo$mag), alpha = 1))

Naozaj to funguje. Tvar súhvezdia, ktorý je bežne zobrazovaný, je zrkadlovo prevrátený (https://upload.wikimedia.org/wikipedia/commons/5/56/Leo_IAU.svg), ale v podstate ide o ten istý obrazec, teda sa nám podarilo vizualizovať jedno súhvezdie. Skúsme to spraviť pre všetky zvolené súhvezdia zo znamení zverokruhu. Ak volíme fazetové zobrazenie, je nutné nastaviť parameter scales ako voľný, pretože inak by sme nevideli nič (zlá mierka).

hygdata_2_vi_z

ggplot(data = hygdata_2_vi_z) +
  geom_point(mapping = aes(x = hygdata_2_vi_z$rarad, y = hygdata_2_vi_z$decrad, size = 1/100^(hygdata_2_vi_z$mag), alpha = 1), colour = "black", background=2) + 
facet_wrap(~con,scales = "free") + theme(panel.background = element_rect(fill = 'white', colour = 'white'))

Ignoring unknown parameters: background

Niektoré súhvezdia sú rozpoznateľné lepšie, iné horšie. Dalo by sa to zlepšiť napríklad iným nastavením filtračného parametra magnitúdy pre každé súhvezdie, podľa počtu hviezd. Niektoré súhvezdia totiž obsahujú aj hviezdy s menšou magnitúdou, ktoré sú ale význačné pre identifikáciu tvaru. My sa však ešte pozrieme na ďalšiu charakteristiku. Pozrieme sa, ako budú tieto súhvezdia vyzerať o 10 000 rokov. Takýto posun môžeme vizualizovať na základe radiálnych rýchlostí za rok, ktoré sú taktiež prítomné v databáze. Takže využime funkcie balíka dplyr:

hygdata_2_vi_z

hygdata_2_vi_z <- (mutate(hygdata_2_vi_z, rarad_10k = 100000*pmrarad + rarad, decrad_10k = 100000*pmdecrad + decrad ))

Teraz už zostáva tvar súhvezdí o 100 000 rokov len vykresliť. (nezabúdame, že je to pri uvažovanej konštantnej rýchlosti pohybu ako aj magnitúde, čo nie je pravidlom):

hygdata_2_vi_z

ggplot(data = hygdata_2_vi_z) +
  geom_point(mapping = aes(x = hygdata_2_vi_z$rarad_10k, y = hygdata_2_vi_z$decrad_10k, size = 1/10^(hygdata_2_vi_z$mag), alpha = 1), colour = "black", background=2) + 
facet_wrap(~con,scales = "free") + theme(panel.background = element_rect(fill = 'white', colour = 'white'))

Ignoring unknown parameters: background

Takto budú naše súhvezdia vyzerať o 100 000 rokov. Tvar niektorých zostane ešte rozpoznateľný, zatiaľ čo u iných bude deformovaný. (Pri parametri 1 000 000 rokov už sú súhvezdia nerozpoznateľné)

V tejto práci sme sa pokúsili využiť programovací jazyk R na organizáciu a vizualizáciu dát. Zamerali sme sa práve na tieto dve oblasti, pretože ide o najčastejšie používané typy práce s dátami. V prvom našom konkrétnom prípade sme formulovali hypotézu o vzťahu magnitúdy a vzdialenosti hviezdy, ktorú sme po preusporiadaní dát a ich vykreslení upravili. Dokázali sme tak, že vizualizácia a organizácia dát je veľmi dôležitá pre sformulovanie/potvrdenie/vyvrátenie vedeckej hypotézy. V druhej skôr zábavnej časti sme sa venovali vykresľovaniu zdanlivej polohy hviezd na oblohe podľa súhvezdí. Zredukovali sme výstupy z databázy, čo nám umožnilo rozumne vizualizovať dáta. Pridanou hodnotou bola prognóza, ktorá snáď budúcim astronómom uľahčí orientáciu na oblohe. :) Metódy by bolo možné zlepšiť použitím komplexnejších postupov.